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INTRODUCTION 


This  report  describes  research  which  was  directed  at  the 
development  of  an  orthotropic  plate  finite  element  for  the 
analysis  of  plate  and  shell-like  structures  which  exhibit 
coupling  between  extension  and  bending.  The  element  is 
especially  useful  in  the  analysis  of  structures  which  are 
fabricated  from  laminated  composite  materials.  The  report  is 
written  so  that  it  would  describe,  for  the  finite  element 
expert,  the  analytical  techniques  utilized  in  the  development 
of  the  element  ana  also  would  be  of  use  to  a  program  "user" 
not  having  the  overall  expertise  of  an  expert. 

The  notation  and  the  methods  used  for  the  definition  of 
element  material  properties  have  been  chosen  as  a  result  of 
a  careful  survey  of  the  literature  on  composite  materials. 

The  notation  and  definitions  chosen  are  considered  to  be  in¬ 
dustry  standard  and  are  best  summarized  in  reference  3,  which 
is  rapidly  becoming  a  standard  text  for  the  analysis  of  com¬ 
posite  materials. 

An  industry  standard  computer  program  SAP  IV  was 

selected  as  host  program  to  accept  the  new  composite  plate 
finite  element.  The  SAP  IV  Finite  Element  Computer  Program 
was  designed  to  easily  accept  new  elements  into  its  element 
library.  The  new  element  rauat  be  self  contained  since  the 
general  philosophy  and  program  structure  is  "overlayed"  into 
the  computer.  The  laminate  composite  plate  element  is  the  new 
element  to  be  integrated  into  the  element  library.  Called 
TYPE  9,  the  new  element  is  similar  to  the  SAP  element  TYPE  6 
in  both  description  and  input.  The  main  difference  is  that 
element  TYPE  9  has  the  ability  to  describe  the  effects  of 
coupling  between  in-plane  extension  and  out-of-plane  bending. 
Element  TYPE  9  is  a  quadrilateral  element  and  is  formulated 
from  quadrilateral  shape  functions  rather  than  from  four 
triangles  as  in  TYPE  6.  Also,  TYPE  9  allows  material  directions 


to  be  arbitrary  for  ease  of  material  input  descriptions.  The 
element  is  modelled  after  the  structure  of  element  TYPE  6; 
therefore  element  TYPE  9  can  degenerate  to  element  TYPE  6. 
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SECTION  II 

LAMINATE  COMPOSITE  FuAT  PLATE  ELEMENT 


The  laminate  composite  fiat  plate  element  is  based  on  thin 
plate  theory  with  the  exclusion  oi  transverse  their  deforma¬ 
tions.  The  following  sections  describe  the  basic  formulation 
of  the  finite  element. 

1.  ELEMENT  POTENTIAL  ENERGY  FUNCTIONAL 

The  principle  of  minimum  potential  energy  furnishes  a 
variational  basis  for  the  direct  rormulation  of  the  element 
stiffness  equations  and  loading  functions.  The  potential 
energy  of  the  element  is  formed  from  the  sum  of  strain  energy 
(U,  )  and  the  potential  of  all  applied  loads  (V  );  i.e., 

?  D 


The  principle  can  be  states  as  follows:  Among  all  the  displace¬ 
ment  functions  of  admissabie  form,  those  that  satisfy  the 
element  equilibrium  conditions  make  the  potential  energy  func¬ 
tional  obtain  a  stationary  value.  Thus, 


6  tt  -  6  U  +  5  V  -  o  ( I ) 

?  P  P 

where  5  is  the  first  variational  operator. 

It  can  be  shown  that 


6U 


p 


i5t  aV 


(3) 


where  o'1  is  a  vector  of  stress  components, 

t  the  corresponding  vector  of  strain  components  and 
V  the  volume  of  the  element. 

Note:  All  vectors  will  be  underscored  with  a  straight  bar  and 

matrices  will  be  underscored  with  a  tilda.  The  super- 
T 

script  of  the  vectors  and  matrices  designates  the 
matrix  is  transposed. 


The  correspond ing  first  variation  of  the  potential  forces 
becomes 

5Vp  =  - J  bT  5u  dV  -  C  £T  6u  dS  (4) 

V  S 


where  b 

vepre 

seats  the 

element 

body  forces , 

is  a 

vector  or 

surface 

tractions  applied  on 

sun  dec  S ,  and 

u 

is  a 

vector  of 

element 

displacements . 

Note  that  the  surface  traction  integral  can  be  used  to  include 
the  point  concentrated  forces  on  the  boundaries  of  the  element. 

The  elements  of  the  strain  potentials  of  equation  (3)  will 
eventually  lead  to  the  element  stiffness  and  initial  load  vec¬ 
tors  and  the  elements  of  the  applied  load  potential  of  equation 
( - )  will  produce  the  various  element  vectors. 

2.  QUADRILATERAL  SHAPE  FUNCTIONS 

The  element  formulation  is  a  geometrically  linear  quadri¬ 
lateral  containing  the  four  corner  nodes  as  shown  in  FIGURE  1. 


FIGURE  1  -  Quadrilateral  Element  Geometry 
a.  Geometric  Shape  Functions 

The  element  shown  in  FIGURE  1  is  described  in  the  local 
coordinate  of  the  element  and  all  material  reference  is  made 
with  respect  to  the  element  local  x  axis.  The  element  area 
domain  can  be  described  by  using  a  polynomial  as 


4 


(5) 


T 

x  =  £  e 

T 

y  =  £  I 

where  x  and  y  are  the  local  coordinates  as  in  the  element 
domain , 


T 

d"  =  Li,  r,  s,  rsj  ( 6  ) 


the  row  vector  of  polynomial  coefficients  jS,  a  vector  of 
generalized  coefficients,  and  r,s,  the  element  natural 
coordinates . 

The  generalized  coefficients  can  be  solved  for  by  evaluating 
the  polynomials  at  the  vertices  of  the  element.  Therefore, 


x 

y 


H1  x 


T 

H  Z 


(7) 


T 

where  H  contains  the  terms  of  the  shape  function  of  the 
element;  x  and  y  are  vectors  containing  the  element  ver- 
rices  s. s  9 


x‘  =  Lx,,  x2,  x  ^  3  x^J  (8) 

T 

y  =  Ly-p  y2>  y3»  y4J 

The  terms  of  the  shape  function  can  be  described  by 

fr  =  ^  (1  +  r^r)  (1  +  s^s)  (9) 


where 


r.T  =  L-l,  1,  1,  -1J  (10) 

siT  =  L-l,  -1,  1,  1J 

define  the  natural  coordinates  ot  the  elements. 
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The 

can 


mapping  of  the  element  geometry  and  displacement  functions 
be  obtained  by  defining  the  Jacobian  transformation  as 


3  =  J  3 

-n  ~  -i 


(11) 


where  J  is  the  Jacobian  matrix  defined  as 


and 


3 

— n 


■*£ 


3x 

1 

sy 


(12) 


(13) 


are  the  first  derivative  operators  in  the  natural  and  local 
reference  frames  respectively. 

Note:  The  " , "  subscript  implies  "partial  differentiation  with 
respect  to".  The  inverse  transformation  is  obtained  by 

=  T*  5  in  (1 


where  J*  is  the  determinant  of  the  Jacobian  matrix  given  as 


J* 


(16) 


and , 


6 


It  will  be  necessary  to  obtain  second  derivates  in  the 
local  reference;  therefore 


4 


E  in 
-  — n 


+  f  al 

-  — n 


(1 


represents  the  second  partial  operator  in  the  local  reference 
given  as 


3x3y 


and  the  natural  set  as 


3 


2 


-n 


/ 


3r3s 


(1 


(2 


The  E  matrix  is  defined  as 


where 


e  .  . 


I  „  T  i  1 

J*  -n  -i 


(22) 


with  gt  being  the  ith  row  partition  out  of  the  G  matrix  and 


8 

-n 


in  -  T*  in  J 


(23) 


The  F  matrix  is  defined  as 


T 

hi 

f  t 
-22 


f  1 
±12  J 


(24) 


where 


with 


f  T 

-ij 


—  2 

m 

li  Ij 

(25) 

gii 

gi2  ° 

(26) 

o 

gil  gi2 

>f  the 

G  matrix. 

b.  In-Plane  Displacement  Shape  Function 

The  plate  element  is  assumed  to  have  in-plane  deformations; 
therefore  the  variation  of  x  and  y  displacements,  u  and  v  res¬ 
pectively,  can  be  expressed  using  the  same  shape  functions  for 
the  geometry  as  in  the  previous  section.  Then 


«T  % 

»T  2v 


(27) 


8 


are  the  domain  displacements  of  the  element 


where 

T 

3U  - 

m 

U1  u2  U3  U4J 

(28) 

and 

3v  = 

U1  v2  V3  V4j 

contain  the  in-plane  model  displacements.  Therefore,  it  is 
assumed  that  in-plane  displacements  vary  linearly  within  the 
element . 

c.  Transverse  Displacement  Shape  Function 
The  plate  element  defined  by  thin  plate  theory  must  have 
a  transverse  shape  function  to  allow  for  proper  bending. 
Therefore  it  is  assumed  that  the  shape  function  polynomial  is 

T  ,  2  232  233  3 

$  =  Ll,r,s,r  ,rs,s  ,r  ,r  s,rs  ,s  ,r  s,rs  J 

(29) 

The  natural  degrees  of  freedom  allowed  per  node  for  bending  are 

2oi  =  Lw  w’r  w’sJ  i  (30) 


Letting 


qQ  =  f  B  (31) 

where  is  the  full  set  of  degrees  of  freedom  in  the 

natural  reference  of  the  element,  ip  is  the  matrix 
$>  evaluated  at  the  nodes  and  defined  in  reference 
2  and  3  is  the  set  of  generalized  nodal  coefficients, 
Then  the  transverse  displacement  w  becomes 

w  =  ^  f1  2o  =  ?T  20  (32) 

where  ip  3  i  r,  the  inverse  of  ip  and 

II  i  : .  Mil-  l  r,.in;vn,:;i'  shape  fu  1 1  r  ■  I  i  <  >  i  ■  ,  both  o  I  which 
,  i  r  ■!  .  h  !  i  1 1 '  •<  I  .n  r  •(  •  I  i  •  r  ■<  •  1 1<  <  ■  7  . 
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The  transverse  displacement  shape  functions  are  defined  in  the 
natural  reference  of  the  element  to  allow  for  ease  of  develop¬ 
ment  since  second  derivatives  must  be  taken.  Since  a  global 
transformation  by  node  is  to  be  performed  at  a  later  stage, 
the  transformation  by  node  from  the  natural  to  local  coordi¬ 
nate  can  be  made. 

3.  PLATE  STRAIN  FUNCTIONS 

The  classical  assumptions  of  linear  thin  plate  theory  are 
made,  essentially  reducing  the  three-dimensional  equations  of 
elacticity  to  a  two-dimensional  set  of  plane  stress  equations, 
for  the  elastic  continuum  of  the  plate,  the  following  assump¬ 
tions  are  made: 

•The  thickness  (h)  is  small  compared  to  the  dimensions  of 
the  plate  in  the  x  and  y  directions. 

•A  line  element  through  the  thickness  remains  normal  to  the 
mid-plane  surface  under  all  states  of  deformation,  inde¬ 
pendent  of  its  translation  or  rotation. 

•The  plate  can  be  isotropic,  orthotropic  or  comprised  of  a 
number  of  orthotropic  laminae,  where  each  lamina  obeys 
Hooke's  law. 

•The  displacements  u,  v,  and  w  in  the  x,  y,  and  z  direc¬ 
tions  respectively,  are  small  when  compared  to  the  plate 
thickness . 

•The  reference  axis  is  taken  as  the  middle  of  the  plate  at 
h/2,  h  being  ’the  total  plate  thickness. 

•The  normal  strain  in  the  z-direction  is  assumed  to  be 
zero,  giving 


e 


z 


w  ■ 


0 ; 


therefore,  the  lateral  deflection  is  given  by. 


w  =  w  ( x ,  y )  . 
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•  St.  Variant's  principle  applies.  That  is,  local  deformation 
occurs  in  the  area  of  applied  loads  while  at  distances 
away  from  the  load,  the  deformation  state  is  not  grossly 
affected . 

•Transverse  shear  deformations  are  neglected; 

Y  =  Y  =  0  . 
xz  yz 

•Displacements  are  linear  such  that 

u  -  uQ(x,y)  -  zw,x  and  v  =  vQ;x,y)  -  zw,  , 

where  w,  =  -6  ,  w,  -  °  ,  and  u,.  and  vn  are  the  in-plane 
’x  y’  y  x  0  0 

displacements  of  the  middle  surface.  The  rotations  about 

the  x  and  y  axes  are  given  by  9  and  9  ,  respectively. 

x  y 

a.  Midplane  Strain  and  Curvatures 

The  mechanical  strains  associated  with  plate  stretching 
nd  bending  can  be  written  as 

e=e+z<  (3j) 

-o 

where  the  mid-plane  strains  are 

(  3h  ) 


the  plate  curvatures  are 

xx  j 

(35) 

yy  ( 

1  xy 


l  i 


with  u  and  v  being  the  in-plane  displacements  and  w,  the  trans 
verse  displacements.  The  thermal  strains  can  be  written  as 


e.  =  a  T  +  zdT 

-t  -  o  -  g 


(36 


where  d  is  a  vector  of  thermal  expansion  coefficients 
relative  to  mid-plane  strains,  Tq  the  element  mean  temp- 


erence  and  the  element  thermal  gradient 


eraturc  dir 
through  the  plate  thickness. 


b.  Strain  Displacement  Functions 

The  connection  between  strain  and  displacement  is  made 
realizing  that  the  in-plane  and  transverse  displacements  have 
been  made  relative  to  a  set  of  nodal  displacements.  Equation 
(33)  can  be  written  as 


£  =  h  Si  +  z!0  30 


(37 


where  8^  is  the  in-plane  strain-displacements  relative  to 
CjT  (the  in-plane  nodal  displacements); 

6r  is  the  transverse  strain-displacements  relative  to 
(the  natural  transverse  nodal  displacements). 

The  in-plane  a isplacements  are 


Si 


(38 


and  the  3-,  becomes 


where 


(40) 


The  out  of  plane  displacements  are 


I  So4  j 

l  j 

where  the  sub-elements  of  the  partition  are  defined  by 
equation  (30)  and  becomes 


(41) 


h  [5  ?n  +  l  ?L] 


(-2) 


where 


f-: 


0 

0  -1 

0  0 


0 

0 

-2 


(43) 


Hx 

~n 


8  R' 

-n  - 


y-T 


(44  ) 


and 


-T 

H 

-nn 


2  — T 

-n  - 


(43) 


Note:  In  equation  (37)  the  z  variable,  which  is  the  plate's 

normal  coordinate,  is  maintained  distinctly  since  it  is 
independent  of  the  in-plane  variables.  Later,  when  the 
strain  energy  is  formed,  the  z  variable  will  integrate 
through  the  thickness  and  merge  into  material  property- 
matrices  . 
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PLANE  STRESS  COMPONENTS 


The  stress  components  for  a  thin  plate  can  be  written  in 
vector  form  as 

o  =  C  (e  -  |T>  (46) 

where  C  ts  one  material  matrix  described  in  the  plate  local 
axes  and  is  expressed  as 

£  =  RT  C  R  (47) 

~ c  "  ~e 

with  C  being  the  material  matrix  in  the  principal  material 
directions  of  the  fibers  and  Ri  being  the  strain  trans¬ 
formation  matrix  from  element  focal  coordinates  to  prin¬ 
cipal  fiber  directions.  The  elements  of  equation  (47) 
are  found  in  reference  2. 

The  dements  of  the  thermal  strain  involving  the  thermal  coef¬ 
ficients  are  defined  as 


a  -  R1  a  . 
~£ 


(48) 


Once  the  material  matrix  is  defined,  the  elements  of  the  material 
elasticity  matrix  can  be  defined  as  in  the  following  section. 

3.  MATERIAL  ELASTICITY  MATRICES 

The  material  coefficients  are  defined  with  the  use  of 
equation  (3)  in  a  slightly  different  form: 


e  dV 


(  h  9 ) 


Since  the  local  z  dimension  is  small  compared  to  the  x  and  v 
plate  dimensions,  it  is  convenient  to  define  the  stress  resul¬ 
tants  and  moment  resultants  as 


A 


adz 


(SO) 


!  4 


ana 


I  - 

J  t 


z  d  z 


rhen,  a  new  stress-strain  matrix  can  be  defined  as 

-  A-  _  f*  ?1  l-o 

2  |m  IbT  pj  (jc 


where 


-JA 

-  fj 

-  I  S 


C  dz 


C  zdz 


C  z  dz 


Letting 


f  A  B 


BT  D 


-T 


equation  (49)  can  be  written  as 


X  2?  5  dv  -  A  |T  1  dA  -  /A  1 


eT  e  jo 

'  -TjTg/dA 

(59) 


L 


where 


AT 

-  A 

I  O  I 

a 

dz 

(60) 

-T 

-  L 

C 

a 

zdz 

2 

(61) 

t- 

Ql 

1  ft 

c 

a 

z  dz  . 

(62) 

The  material  A,  B,  D  matrices  and  the  thermal  load  coef¬ 
ficients  A^,  B  and  D  can  be  related  to  laminar  material  by  posi- 
-  *  -T  “1 

tion  t  in  the  material  build-up  as 


A  = 

L 

z 

i  =  l 

c.  (t.-t.  ,) 

~i  l  l-l 

(63) 

B  = 

1/2 

L  0 

&  ^  ^  - 

(64) 

D  = 

1/3 

L  -  5 

k  k  ' 

(65) 

A„  = 

- 1 

L 

z 

i  =  l 

Si  Si  (ti  -  'i 

(66) 

-T 

1/2 

L  ~  ~  , 

Y  C-  a.  (tf 

1 

-  'Id 

(67) 

-T 

1/3 

y  C,  a-  (t? 
i  =  l  1 

-  *Lx> 

(68) 

where  the  subscript  "i"  implies  coefficient  evaluation  at 
laminae  level  "i"  and  L  is  the  total  number  of  fiber  lamina 
levels . 

6.  COORDINATE  TRANSFORMATIONS 

The  element  information  is  initially  determined  in  the 
natural  coordinates  of  the  plate  since  it  is  quite  easy  to  expres 
all  loading  and  stiffness  information  in  that  reference.  Ulti¬ 
mately  the  information  must  be  transformed  to  local  coordinates 
(x,y,z)  and  also  to  global  coordinates  (X,Y,Z).  The  following 
sections  describe  the  transformations. 


•-i.  .Natural  to  Loral  Transformation 

The  natural  coordinate  variables  per  node  are  defined  as 

/  \ 


j  *  1 

\  W’r  /  i 

!w’s 

i  I 

where  6  is  the  rotational  degree  of  freedom  normal  to 
plate  at  node  "i" . 

he  transformation  matrix  required  becomes 


tn(i,i  ^Sli 


where 


t-  .  =  0 

~nKi 


I  m 

L  O1  0:  1_ 


L  u  V  w  u  G  0J 
x  y  z 


with 


i:i!  .■  *  i 


.Noting  that  the  rotation  degrees  of  freedom  are  definec  as 


-W’x  I 


the  v  matrix  becomes 


Q; 


(74) 


The  complete  natural  to  local  transf ormation  therefore  becomes 


~  n£l 


~n£  2 


-n£  J 


:n£4 


(75) 


b.  Local  to  Global  Transformations 

The  local  to  global  transformation  quantities  are  somewhat 
more  difficult  to  obtain  since  the  transformation  involves  the 
local  coordinates  of  a  quadrilateral  element.  Obviously  only  3 
points  define  a  plane;  the  fourth  point  of  the  quadrilateral  is 
unnecessary.  The  fourth  point,  however,  may  not  lie  in  the  same 
plane  as  the  other  three  points.  Therefore  local  transformations 
by  node  are  determined  and  are  averaged  to  obtain  a  general 
transformation  used  in  determining  the  element  coordinates  and 
m  transforming  element  matrices  where  applicable. 

Defining  the  nodes  of  Figure  1  as  i,  j,  k  and  1  and  allowing 
this  sequence  to  permute,  the  element  normal  coordinate  at  node 
"i"  which  also  permutes,  is 


zi 


V  .  .  X  V  - , 

Zli _ Til 

IV  .  .  X  V  .  J 
-li  -il 


(76) 


where  the  "X"  symbols  denote  a  cross  product  of  two  vectors, 

and  V . ■  implies 
-li 
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I 


z.  -  z.\ 

3  1) 


v. i t r.  X,  Y  and  Z  being  global  coordinates. 

A  global  material  reference  is  defined  as  a  global  vector 
defined  as  input  which  locates  the  general  local  "x"  of  ail 
elements  referred  to  that  vector.  All  material  properties  are 
defined  relative  to  this  x  coordinate  which  becomes  the  ele¬ 
ment's  local  x  coordinate.  The  element's  local  y  coordinate 
at  each  node  can  be  calculated  as 


e  .  X  V 
zi  -m 

e  ■  X  V-] 
zi  -m 


(7  c) 


’inally,  the  local  x  coordinate  at  each  node  is  determined  as 

(79) 


§  .  X  g  . 

yi  21 


XI 


j€  ■  X  §  . 
i  yi  zi' 


The  local  to  global  transformation  at  node  "i"  becomes 


reT  "1 
x 


„  i 

ey 


(SO 


L&: 


where  the  direction  cosines  of  each  coordinate  are  placed 
in  row  order  in  the  transformation  matrix. 

The  average  transformation  used  to  determine  the  local 
coordinates  is  obtained  by  first  averaging  the  nodal  normals  as 


e 


z 


(81) 


and  substituting  into  equations  (78,79  and  80)  to  produce  an 
average  t£g. 
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"he  global  degrees  of  freedom  can  be  defined  by  node  as 

'u  ' 

V 

■gi  <82> 

e, 

8, 


where  U,  V  and  W  correspond  to  global  displacements  rela¬ 
tive  to  X,  V  and  Z  respectively  and  0‘s  corresponds  to 
global  rotations  about  X,  Y  and  Z  respectively. 

The  complete  transformation  becomes 


where 


~lgl 


:ig2 


~<Lg  3 


’  g4 


(83) 


.igx 


(84) 


and  q „  and  q  are  the  complete  list  of  degrees  of  freedom 

x-  -*g 

per  element. 

If  the  quadrilateral  element  is  perfectly  flat  in  its  space, 
then  the  t.  •  becomes  exactly  t.  .  If  it  is  not,  then  the 

-£gi  7  ~«,g 

element  space  appears  as  a  curved  space.  This  effect  should 
allow  the  element  to  behave  as  a  shallow  shell. 


c.  Natural  to  Global  Transformation 

The  element  stiffness  will  be  transformed  from  natural  to 
global  coordinates  directly.  Therefore  that  transformation 
becomes 
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£ 


n 


(85) 


where 


~ng 


(86) 


Since  both  original  transformation  matrices  are  partitioned 
diagonally,  the  nodal  transformation  matrix  is 


t  .  =  t  „  .  t .  . 

~ngi  ~n£i  ~£gi 


(87) 


which  contains  many  off-diagonal  zeroes.  Therefore  it  is  con¬ 
venient  to  define  transformations  by  displacement  and  rotation 
degrees  of  freedom.  That  is,  the  displacement  transformation 
at  node  i  is 


!Di 


!i 


(88) 


and  the  rotation  transformation  at  node  i  is 


This  modification  will  save  transformation  operations  later  on. 

7.  ELEMENT  STIFFNESS  MATRIX 

The  element  stiffness  is  easily  defined  in  the  natural  coor¬ 
dinates  of  the  plate,  given  the  strain  displacement  functions. 
Once  this  stiffness  is  determined,  it  can  be  augmented  with  a 
scaffolding  or  artificial  torsional  stiffness  for  the  plate's 
normal  degrees  of  freedom.  Finally,  this  stiffness  can  be 
transformed  to  local  and  then  to  global  coordinates  for  assembly 
into  a  master  stiffness  matrix.  The  following  details  the  above. 
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a.  Plate  Element  Stiffness  in  Natural  Coordinates 
Defining  a  new  set  of  degrees  of  freedom  as 

5-  ■& 

the  strain  components,  as  defined  by  equation  (58)  become 


B 


~n 


9 


n 


(91) 


where 


(92) 


The  second  integral  of  equation  (59)  can  be  used  to  define  the 
element  stiffness  in  local  coordinates  as 


E  e  dA 
~m  - 


K  q 
~n  -n 


where 


(93) 


R 


~n 


E 


~m 


B  dA 
~n 


(94) 


For  convenience  or  computation,  the  material  matrix  Em  in  equa¬ 
tion  (94)  is  Cholesky  factored  as 


T 

E  =  U  U 
~m  ~ 


(95) 


where  U  is  a  upper  tri-diagonal  factoring  matrix. 

This  allows  equation  (94)  to  be  written  in  a  more  efficient 
form  as 

K  -  f  (U  B  )T  (U  B  )  dA  (96) 

~n  J  ~  ~n  ~  ~n 

A 

which  allows  the  triple  matrix  product  to  be  replaced  by  a  sim¬ 
pler  transpose  symmetric  product.  This  process  is  especially 


2? 


efficient  since  the  Choi « sky  factoring  is  performed  (at  most) 
once  per  element.  Also,  numerical  integration  is  to  be  per¬ 
formed.  Savings  will  occur  at  each  integration  point  after 
the  first. 

b.  Artificial  Torsional  Stiffness 

The  flat  plate  theory  does  not  have  any  mechanism  to  di¬ 
rectly  include  twisting  of  the  plate  normal  to  the  plate  sur¬ 
face.  Therefore,  if  two  coplanar  elements  are  assembled  at  a 
common  node,  a  singular  stiffness  exists.  To  avoid  this,  an 
artificial  or  scaffolding  stiffness  is  added  to  the  normal 
rotational  degree  of  freedom  6  .  There  is  no  change  in  the 
system  equilibrium.  For  convenience,  this  is  performed  at  all 
nodes  of  the  element  since  it  would  be  difficult  to  determine 
coplanar  effects  in  general.  This  does  change  the  overall 
element  equilibrium.  If  the  amount  of  artificial  stiffness  is 
xept  small  and  the  local  rotational  stiffness  effects  are  in 
equilibrium,  then  the  error  can  be  minimized.  Defining  a  vec¬ 
tor  of  normal  rotations  at  the  nodes  as 


the  artificial  torsional 


'fness  matrix  relative 


to  tj  becomes 


where  f  is  an  input  scaling  factor  which  can  vary  as 


0  <  f  <  1 


(99) 


and  can  be  set  in  l.E-8  increments,  C  is  an  artificial 
coefficient  estimated  from  element  bending  stiffness 
coefficients  and  element  area;  i.e., 

C  =  MIN  (D( 1 , 1 ) ,  D( 2 , 2 ) )  *  AREA  (100) 

with  the  D's  defined  in  equation  (55). 
c.  Natural  Stiffness  Matrix 

The  degrees  of  freedom  q  and  6^  defined  by  equations  (90) 
and  (97)  can  be  merged  to  degrees  of  freedom  described  in 
equation  (69).  This  requires  the  re-ordering  of  stiffness  coef¬ 
ficients  of  equations  (96)  and  (98)  to  produce  a  natural  stiff¬ 
ness  matrix: 


K 


~n 


MERGING 

REORDERING 


(101) 


relative  to  q  . 

-»n 

d.  Global  Stiffness  Matrix 

The  global  stiffness  matrix  K  is  formed  by  transforming 

o  ^ 

K  from  natural  coordinates  to  global  using  equation  (85).  The 
transformation  is  formed  using  equation  (1)  realizing  that  the 
strain  energy  U  is  invariant  relative  to  any  coordinate 
reference.  Therefore, 


U 


P 


Using  equation  (85)  produces 


K 

~P 


-HP. 


(102) 
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Tin  triple  matrix  j  -  r>  •<  tin  I  i  m  |  >  I  i  <  -i  i  in  <  <|  u.l  I  i  on  (  I  0  t  )  i  :  quite 
inef  f  icicnt ,  especially  since  T  %:■  highly  d  i  agonal  .  Eft  .  c  i  cncy 
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can  nowtv.T  s>e  eincua  sy  partitioning  k  and  K  into  i  X 

~  n  ~  p 

n  ?  & 

sub-.natr  ices  laoeliea  K.  ..  ana  K„.  y*  where  i  and  i  range  from 


to 


?hen 


-  J 


t  ■  K  .  ■  t  . 
~l  .13  ~3 


(104  ) 


wnere  t.  matrix  relates  to  equation  (33)  when  i  equal 

1,  3,  5,  7  ana  relates  to  equation  (89)  when  i  equals 

2,  4,  6,  3. 

Additional  et  t  icioncy  is  obtair.ee  wnen  tr.e  triple  matrix  prcuuc^ 
is  performed  such  that  the  right  portion  matrix  multiply  is 
first  formed  ana  positioned  back  into  K.?.  Then,  the-  left 

^  j 

multiply  is  formed  with  the  resulting  product  ana  placed  into 


?.  numerical  integration  of  area  functions 

The  elements  in  the  matrix  of  equations  (94)  and  (9b)  are 
very  difficult  to  integrate  exactly,  therefore  approximate 
numerical  integration  can  be  performed  with  sufficient  accuracy 
for  j  ustif ication.  Gauss- Legendre  Numerical  Quadrature  has 
seen  selected  to  perform  the  integration  of  the  stiffness  coef¬ 
ficients  as  well  as  other  area  functions.  Equation  (94 )  can  be 
rewritten  and  transformed  relative  to  variables  and  limits  of 
..  *  Z  Ml  £  r  itlJuo  of  S 


I  f.  BT  (r,s)  E  r  (r,s)  J*  era 
-x  ~n  ~m~n 


(  -  .  i  ) 


An  element  of  this  matrix  c.m  be  written  as 

( 1  f1 

k  .  .  =  J ,  J .  :  .  •  ( r , s )  ends 

l  :  -  i  - 1  1 3 


( I  G  c  ) 


where 


f  ,.•(:■,  s )  =  E  b;,'  (k,i)  I  r(K,n  :-■(«.,:,) 


~n 


(.07) 


I 


The  stiftness  coefficient  can  then  be  approximated  as 

--i£  %  hj  <vv  \  w«  (i 

where  n. ,  n„  are  the  number  of  Legrendre  root  evaluation 

X  4 

points  in  the  r,s  directions  respectively. 

r..  and  s .  are  the  roots  of  the  Legrendre  polynomial, 

Jc  X. 

and  are  the  appropriate  Gauss  weighting  factors. 


or 


t  or 


rk  =  sk  '-yJT  n/T 
Wk  =1,1 


(109) 


n  =  3 , 


Wk 


vi~  vr 

5  ,  0,  5 

5  8  5 

9  ,  9,  9 


(110) 


9.  ELEMENT  MASS  MATRIX 

A  consistent  mass  matrix  relative  to  the  in-plane  variables 
in  one  coordinate  can  be  written  as 


L 


p  t  H  H  dA 


(111) 


and  a  corresponding  lumped  mass  matrix  can  be  formed  by  summing 
the  rows  of  the  consistent  mass  matrix  as 


M 


i 


p  t  f  HT  dA 


assuming  t  and  p  constant  and  realizing 


(112) 


2  G 


The  c0ri7.cner.tL.  o  M"  ire  applied  to  the  translatory  degree..  of 
-•  reedcr.  in  ill  j.ocu  1  dire  cions ,  per  node.  A  rotary  inertia 
ciieet  can  be  included  by  an  approximation ; 


me  gioDax  camped  mass  veccor  can  ne  formec: 


(113) 


a 


1  ;*kp  o  x 

l  :  ‘ 


(114) 


where  a  is  a  o  ±  unit  vector. 

1C.  ELEMENT  LOAD  VECTORS 

The  element  load  vectors  are  established  from  element  prop¬ 
erties  such  as  material  constants,  temperatures,  pressure,  mass, 
area  and  acceleration  constants.  The  following  sections  des¬ 
cribe  the  load  vectors  ‘developed. 

a.  Tkerme 1  Loud  Vector 

The  thi.ro  integral  of  equation  (53)  can  be  usee  tv  def  .ra¬ 
the  element  t normal  load  vector  using  equation  (SI): 


fv 


t  (dA  =  in 


5 “  dA 
~n 


(M 

,j  (xg) 


There f ore , 


(liC) 


is  the  therm.il  vector  relative  to  the  natural  coordinates 


a,.,  i:  I  h. 


I  .  Vi 


mi  ■■  I  1  r  i :  l  ruin  .  : 


i  7 


(II?) 


§t  =  -  k!  ■ 

The  global  thermal  vector  can  be  determine^  by  applying  the 
natural  to  global  transformation. 


T  "  F" 
~ng  -n 


( 1 1  o  ) 


b.  Pressure  Load  Vector 

The  pressure  load,  by  node,  is  formed  for  the  shape  terms 
associated  only  with  the  geometry  and  is  appliec  in  the  2  (plate 
normal)  direction.  Therefore, 


f  p  K  dA 

Ja 


(115) 


A  pressure  load  vector  can  be  formed  relative  to  local  coordi¬ 
nates  as 


A? 


(120) 


The  global  pressure  loan  can  be  formed  as 


^  a  FP- 

■  £g  -  £ 


(121 ) 


which  transforms  the  normal  traction  into  a  global  traction, 
c.  Constant  Acceleration  Loaa  Vector 

The  acceleration  vector  can  bo  computed  from  the  mars,  w-cloi 
defined  in  equation  (114);  i.e.. 


(  t  2  ?  ) 


section  : r i 

MODI:'  ICATT ONS  TO  SAP  IV  COMPUTER  PROGRAM 

The  SAP  TV  of  mputer  program  modified  to  accept  the  new 

composite  element.  The  following  sections  describe  changes  to 
the  existing  progr-am  as  well  as  new  routines;  of  the  composite 
[-.left  e  Lemon  t  . 

1.  SAP  IV  STRUCTURE 

The  main  changes  to  the  SAP  IV  program  occur  in  the  el erne n 
library  control.  They  occur  in  the  element  generation  portion 
of  the  program.  Figures  2  and  3  depict  the  major  routines  used 
:  y  SAP'  to  perform  a  static  analysis  of  its  elements.  The  EI.TYP 
i  .i  stine  calls  a  new  routine  called  PLATE  ,  which  is  usee  tc  ci.  . 
(VEF.LAY  (8).  (stress  recovery)  Once  OVERLAY  (8)  is  callec,  a 
call  is;  made  to  CPLATE,  the  main  routine  of  the  composite  plate 
finite  element.  Figures  4,  5,  6,  and  7  describe  the  flow  of 
routines  used  by  CPLATE. 


•  l 


I 


description 


Input  title  ana  master  control  rare 
Input  nodes 

Write  8  ID,  X,  Y,  Z,  1 

Call  element  type  (new  composite 
plate  element  called  here) 

WRITE  1  Element  stress  matrices 
for  recovery 

WRITE  2  Element  stiffness  matrices, 
mass  and  load  vector:., 

Input  Concentrated  i.  oads/masceu 
WRITE  3  Concentrated  load:-;/ 
masses  for  assembly 

READ  2  Read  stiffness  ror  asoemblv 
WRITE  4  Assembled  Stiffness  and 
Loads  in  blocks 
WRITE  9  Assembled  Mass  Vector 

Call  appropriate  solution  type 
0  Static 

1  Nodal  Extraction 

2  Force  Response  Analysis 

3  Response  Spectrum  Analysis 

4  Direct  Integration 


BASIC  PROGRAM  PLOW 


U.r. 

Static  v.oiut  ion  option 

Solution  or  displacement  equation 
READ  4  Read  Stiffness  and  -o ado 
for  decomposition 
WRITE  2  Solve  for  displacfemer.tr. 

READ  a  Read  ID  identir icat ion 
matrix 

READ  2  Recover  displacements 
Print  displacements 

READ  1  Element  stress  matrices 
Stress  recovery  and  Printing 


C  SOLUTION 


AND  RECOVERY  STAGE 


DESCRIPTION 


Composite  plate  routine 

Input  material  Properties  in 
tables 

Input  global  material  vectors 

Input  load  factors  TLO  as  a 
unction  of  thermal.  Pressure 
and  constant  acceleration 

Input  element  definition 

Form  ABD  matrices  and  calcu¬ 
late  thermal  stress  vector. 
Kodify'E  =  US'  U 

Determine  local  to  global  trans¬ 
formation  and  locax  element 
coordinates . 

PROSN  -  forms  vector  cross  product 
SI’RANS  -  special  transpose 
KG  I.  -  used  to  transform  global  to 
local  coordinates 

l 

Set  element  K  ,  B  ,  A ,  N  matrices 
_ _  ~ n  -n’  -g 


FIGURE  4  FLOW  DIAGRAM  TOR 
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Loop  on  NC  (Gauss  points)  in 
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*  (ii,  :?) 

~n 

Accumulate  B  over  all  Gauss 
-  ,  ~  n 

points 

7KULT  -  forms  ..pec la  1  multiply 
of  L;  h,- 

ATA  -  form  special  symmetric 
muitir  ly  Co  )  •  ( U  •- n ) 


i  >. ) I  'iTi i_>  i-Ti  i ■  _i. 

mat 

rix 

prcauct  «nc 

stores  result 

i  n 

K 

^  i'. 

Lower  symmetr 
forma  K  (Ii, 

12) 

U  C  v : 

^  v.  *  r  ■  it  u  1 1  i  p  i  */ 

W(I1)*W(I2)*K 

Lower  symmetr 

ic 

auc 

-  acs  inte- 

oration  resul 

L  U 

to 

accumulatv r 

for-  local  clement  stiffness. 

Let ermines  area  distribution  vec¬ 
tor  A  by  node . 


1  lAGRAK  FOR  GFLA7F  (i  I'T) 


liiXJ 


.-  pec :..i .i  matr.x  mu:;  ipiv 

u,  -  u 

bn  ~n  ‘-i! 

‘form  element  .  trees  mat 
rolat  Lve  to  .  nii.jV 

'..0 l  location  v e  .toe  oi 
.’•umbers  for  element  pi o 
connectivity . 

•.a  iCulato  banc  w  i  c t n  c  > i 
st ii i ncib ,  x o a 
an n  mans  vecto 
assembly 

stress  matrix 
recovery 

c on t rob  t o  mu i 


'  I  Ax.  RAM  FOR  Cl  LAV:  (CbVi'V) 


WR  I.  T  V  1 


WRITE  1 


I 


Its.-  a  dettori  ption  c.  the  routines  used 
. _tt  riute  overlay. 

. 1  iue  ■  .ilea  by  C PLATE  to  fill  tr.e  A,  b  and  D  o: 
S3,  :  and  55  material  Hid  trices  :  rom  the  input 

•  -  ’  'em:;  arrays.  Equal  ions  CO,  Cl,  and  G2  are 

.  .  :  !  ■  rh*  thermal  load  vectors,  and  thermal 

:  .  .  equation  i; 7 . 

a.i.ae  cdileu  by  L'KA  7  to  idcio:  die  material 
into  an  upper  triangular  matrix  ,u  described  in 


-  .  r  ubrautin.  called  by  CP LATH  to  determine  the  local  to 

•  1  :  carnation  matrix  as  described  by  equations  76 

.  ..  .  line  routine  checks  for  proper  area  definitions, 

.  -  :  ..  n  urine  called  t»y  LOCAL  to  perform  vector  cross 

:  .  .  c :.  ..  in  equations  76,  78  and  79.  The  resulting 

.  .’ir  component,  are  normal  i  nod  to  unit  vectors, 
i  -  ..  .’uji  ouc  in"  Called  by  LOCAL  to  perform  an  in-place 
can.-.:  e  matrix  transpose  of  the  local  to  global  transfor¬ 
ms  lot.  ;•  eq..a  t  i  on  «0  . 

-  -i  u:  :  u'dx  ..led  by  LOCAL  to  transform  global  element 

:•  :  it  i  ■  Icon  element  coordinates .  Trans  format  1  on 

tin  lx  i  u':.\  by  equation  60. 

-  ..  s  u  t-ro ,  i  ?.  Luc  oalle:  Tv  CPI. ATI!  to  .ct  matrix  array  space  to 


a  r  1  \  v  j  .  u  a  . 


:  i cal  ly,  it  is  a  red  to  set  matrix  space 


.■  -  soutine  called  by  CPL.ATE  to  :  orm  composite  plat « 

.  i. train  uisnla  "onent  :  unction  ••  us  described-  hv 

- 1 1 

■  _u.it.  :  ■  del  anon  equ.it  t  ons  ''  t  L ■  -ugh  5  . 

i  -  a  :  s.  it  in-  on  Hud  bv  CLuATf  t,.  ,u  c  ijr.-ilate  1  :.*•  r-  .i,n 

* . ;  n  am  r  .  ...  o  ,  ..  iadraturr-  ...;  uA  --ver  al  :  (juus...  *'■  .ut.  C 

the  are  i  •.  ;  in  equation  ill. 

-  .  ,  ..  ;  '-nil.  !  Is  :  1  ATT.  t  .-  :  or::,  a  :  ,  ■  •  rr. 

the  a..'  ■  i  :  ;  r  oc.ot  ler.erd  bed  it.  eq  utisn  9*4  as:  ■' 

'I’ri  i  .  :  .  no.  I  gin  1  t  o r  ■ ' :  !  i  o  1  s o ■  . ; : , r  1 1 s'  ■*  s •  1  n  t .  - 

.•rat  i  •  ;  .f.  .i  •  ■  ; u r-  • . 


.  nr- 


■  •  -  ijii'y  a*.  .  cm. i  itO  jy  oUauRN  "to  peri orm  a  spec  iai 

■■  i t  r  : x  mu .  r  ;  p  1 7  where  the  leading  matrix  is  an  upper 
Tri..;i.;tuar  .aerix;  described  ay  equation  v6 . 

1  Mibrour  me  oal  2  oh  by  UbrCEKL  to  perform  a  special  svm- 
uc  ::;ul  tipiy  aeeooa  in  equation  96. 

■  ■  -  •  utibro,.  •  stc  culleu  by  C PLATE  to  perform  a  lower  sym¬ 
metric  ::.u  1  tipiy  o:  the  natural  stiffness  components 

-.urine  •  .crriiUi  integration,  as  in  equation  108. 

-  i  s..!  r-.-.s  cabled  by  C PLATE  to  per r orir.  the  lower 

'vj:t:.e- .-t-  - :  r  1  or:  o:  tne  stiffness  matrix  components 

. .1  integration ,  as  described  by  equation  ICS 

-  a  an:  .a..:  culled  by  C PLATE  to  determine  the  distri¬ 

buted  c  ;  a : e a ,  by  noce,  for  the  quadrilateral  plate  elu¬ 
te  n „  .  ini  area  function  is  needed  in  equations  112  and  1 

-  a  s. ,brou: ine  called  by  CP LATE  to  form  a  transformation 
tat: a  x  fr  r..-.*ura_  to  local  coordinates  as  described  in 


:unr it  ir.e  cal  iced  by  CP  LATE  to  form  the  natural  to 
global  coord  !  r.cite  1  runs  format  ion  as  shown  in  equation  67. 
hi  roe  t.  o  t on n formation  matrix  is  diagonal,  only  the 
a idrccidi  are  stored. 

■  -  a  ..  rou :  s.c  .  .uicc  by  CT  LAT..  to  ore:  a  net  :r.i . 

'.lo.  a.  t :  .ns.  rm-s,  i  :r.  of  the  .  tinin-u  i  s;  i.  acer.er.t  .eid  .  :  i\-.  • 

nic.  .  ..c  r.i  nuij’ito.-,  scribed  in  equations  116  ana  12  6. 

.  sn.  c  iiwc  b.y  ■  .  -ATE ,  EOaP.T  ,»nu  :  d.-'.Xbii  to  .  .  :  - 

1  ore  c-.  ge i .c :  MHUX  mu  a  t -.-p..  teat  ton  01  .-irn  .t Sea rv  rs.s  ni  . 

-•T'  Ctfd  :  ret  or  por.  at  ion  to  .any  sub-matrix  position. 

—  .  s  ns'os  t  i.i.e  caiieU  ny  f  fl./.r.1.  to  caaou  1  ate  tr.o  caes.*  ‘  : 

sue...  j.jiv  ,  ccr.tr.inr  acceleration  and  mass  vectors  desert;, 
i  y  a  ,uat  i  119  a  no  122. 

-  1  ...ibrout  j  c.!..es  by  C  F  LATE  t--.  prepare  the  elements  .  : 
*.!•■  .  .•■•me:.;,  p.oral  stli  ■  ne:s  matrix  with  element.,  iron  r.« 

s  j(  .4  i  Si  ..no  art  1  :  sial  tors  i>  siu  a  hit:  :  tness  . a  — 

t-  r  i  -e  ,.  In.  .  ..»r  r>.  ■  •  s  s  is  .K  oOfiLcd  :  r.  •  at  in:,  .is. 


•-  or  the  natural  sr;r: nest  matrix  into 
'loL-al  stif  t  ne:...-,  matrices.  The  original 
■a  matrix  ana  the  final  global  sti fines 
::.e  came  matrix  space.  This  procedure  i 
’nation  in 3. 

r.illed  by  I'TRANbr  to  perform  a  special 
irodacr  used  in  sti !  f  ness  trans  forma  t  ion 
:  r  forms  an  el:  Lcient  doul  -  matrix  mo  it  1 
■*ite  of  t:.e  original  sub-matrix  as  desc; 


called  bv  PTRANSF  to  move  the  elements 
any  rank  and  place  them  into  new  matrix 


called  by  f".  RANSF  to  move  elements  sir. 
that,  the  receiving  sub-matrix  is  the  m 


me  called  bv  CP LATE  to  perform  a  specie 
!  such  that  the  post  multiplying  matrix 
;  shown  in  equation  126. 


ON  IV 

R  i;' I  CATION 

■  s..  s v**'  the  s-s  *  .  program  proOuC'  a  new  ve r'- 
ieh  bar  I  •  iilv  includes  an  additional  element 
.  i :  v.  ukuient  >  'I  he  r.ri  ginai  SAT  IV  program, 
..  r  e  :  v  moat  a  ana  many  di  fretvnt  static  and 
; . .t 1 vs i r  .  .  ocedures ,  has  been  a  ±  r  e  a  d  y  v e r i  — 

.■  the  new  element  was  check c  d  under  various 
s  conditions. 

t  was  compared  with  many  simple  degenerate 
sblems  ana  was  found  to  produce  excellent 
rest ..  were  per  formes,  for  element  TYPE  6  and 
were  more  favorable  for  the  simple  cases. 
■.Sait.:  in  larger  displacements  while  TYPE  6 
.ewer  than  exact  solution  values. 

••ectior.  describes  a  group  of  verification 
.  t  voc  s n e ..  j-  anc  a  doubly  curv'6c  blade.  ....  a 
obtained  bv  classical  plate  and  shell 
: oat  tc  material  properties.  The  pro.  eh  art 
if)'  series,  so  approx.mnte  to.  ut  i  or.s  of  vano 
v  coi.ait  sons  .  once  t n»-  A,  o ,  a  natr .  c.es,  v.ss 

*.""•/  C  ’.Ai'ii  t  1  i.-i  i  S  cij-vi  itit.  f  '.i*-  .j  Yj  j.-T'O  a  -  i.4t;  wit  u 

ci  v t t .  c  v  *;  i  r. h t  - 1  rr.i  vxp  .  Mcs T 

c. :  i  . '  «.•  lornorit ;'i  .-cr  r. j  i\  ct  ion.  rcr.vt.-r^  c.  r:c  <.• 

"  .-.t  .nice  s  see  g.c.eri.  ..teev  tr.c  O-a...- 
;  erer.ee  I. 


4.  L02/902Jt  UNDER  UNIFORMLY  DISTRIBUTED  LATERAL  LOAD 

WITH  ALL  EDGES  SIMPLY-SUPPORTED 

(i)  Size  of  laminated  plate: 

a  =  10"  b  =  10"  t  =  4h  =  .2" 


(2)  Properties  of  plate: 
p  =  .0275lb/in3 


.  0275 
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(3)  Loading  condition: 

p(*,y>  =  PQ 

(4)  Boundary  conditions: 


x  =  0  ,  a ; 

w  =  0 
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FIGURE  11 


[  0 / 9  0  3 1  UNDER  FREE  VIBRATION  WITH  ALL  EDGES  SIMPLY-SUPPORTED 


( 1 )  Size  of  laminated  plate: 

a  =  10"  b  -  iO"  t  =  2h  =  .2" 

(2)  Properties  of  plate: 
p  =  .0275  lb/in3 
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(3) 

Loading  c 

ondition : 

free  vibration 

(4) 

Boundary 

conditions : 

lb 


x  =  0,  a; 
y  =  0,  b; 

requency : 
(f ) 


CXciC  l 


(f ) 


6.990  Hz . 
8.867  Hz. 


6w  =  0 , 
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[SO/ 0 9/90]-  CURVED  PLATE  UNDER  UNIFORM  PRESSURE  WITH  ALL 
EDGES" SIMPLY- SUP PORTED 


tn  w 


(1)  Type  of  Blade:  (Figure  16) 

A  typical  blade  used  in  aircraft  engines,  named  J79  B/AL, 
was  approximated  using  64  quadrilateral  plate  elements  in  SAP4A 
(both  TYPE  6  and  TYPE  9).  An  equivalent  NASTRAN  model  was  also 
run  along  with  an  experimental  test  to  determine  fundamental 
frequency  at  zero  frequence  (results  from  AFSC  -  Wright-Patterson 
AFB) . 

(2)  Properties  of  blade: 

Material  used  corresponded  to  the  input  used  in  the  NASTRAN 
run  using  anisotropic  material  properties, 
leading  edge: 

;  =  .000407  lb  sec in  . 

Cxx  =  26.9E8  Ib/in2  Cxy  =  4.6E7  lb/in^ 

Cyy  =  20.3E8  lb/in2  Gxy  =  7.3E7  lb/in2 

Blade: 

p  =  .000251  lb  sec^/in 

jd me  material  coefficients  as  above. 

(3)  Loading  Condition: 

Mass  and  Stiffness  distributions  for  eigensolution . 

(4)  Boundary  Conditions: 

Base  of  blade  completely  fixed  and  rest  of  blade  free. 

(5)  Natural  Frequency  of  Blade:  (Fiest  FLEX) 

(f)exp  =  110  Hz. 

(  i  ) Nar.tr  an  -  ;  0  6  .7  Hz  . 

( f )sap  (T6 )  =  i 08 . 3  Hz. 

( f ) sap  (T9)  -  L 1 2 . 8  Hz. 
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SECTION  V 

DISCUSSION  AND  CONCLUSIONS 


A  quadri lateral  composite  plate  finite  element  has  been 
a -idea  to  the  SAP  IV  computer  program  library  to  be  used  on  plate 
t'.'pe  structures.  The  element  is  a  four-noded  sub-parametric 
; lat  plate  element  excluding  transverse  shear  deformations.  The 
element  is  "incompatible"  relative  to  mid-plane  surface  rota¬ 
tions  along  the  inter-element  boundaries.  The  element  converges 
relatively  well  despite  the  incompatibility,  as  long  as  the 
clement  maintains  a  relatively  rectangular  shape  with  a  reasonable 
element  aspect  ratio. 

The  element  was  run,  modelling  various  simple  plate/beam 
configurations  and  showed  excellent  results.  More  complex  models 
were  designed  to  test  the  composite  laminate  behavior  of  the 
.foment,  as  described  in  section  IV.  The  models  included  flat, 
cylindrical  and  doubly-curved  shells.  The  results  obtained 
compared  favorably  with  classical  series  solutions,  (approximately 
1-6%  disagreement)  The  results  obtained  in  section  4.9  for  the 
typical  blade  conf iguration  shows  that  the  NASTRAN  and  SAP 
(TYFE  6)  element  values  were  slightly  lower  than  the  experimental 
number  while  the  result  of  the  new  element  (TYPE  9)  was  sligntly 
higher.  The  new  element  is  based  on  an  incompatible  formulation 
and  in  general,  convergence  is  not  guaranteed.  But,  in  more 
erases,  the  element  has  been  found  to  be  slightly  stiffer  than 
compatible  elements  and  therefore  the  frequency  is  higher. 

The  quadrilateral  element  will  be  inherently  stiff  if  the 
four  points  do  not  represent  a  flat  surface.  Therefore  the  "quae" 
element  was  relaxed  to  better  represent  shell  behavior  by  allowing 
the  flat  plate  element  stiffness  coefficients  to  be  transformed 
relative  to  the  local  nodal  coordinates  oi  the  element. 

The  thin  plate  theory  used  to  develop  the  element  does  not 
account  lor  normal  torsional  effects.  There 1  ore ,  non  global 
adjacent  elements,  when  assembled,  will  produce  a  singularity 


normal  to  the  plate  it  the  two  elements  are  coplanar.  To  avoid 
this  singularity  in  general,  an  artificial  torsional  stiffness 
or  scaffolding  matrix  was  added  to  each  of  the  four  nodes. 

This  does  violate  element  equilibrium  but,  if  the  magnitude  of 
toe  coefficients  are  maintained  relatively  "soft"  compared  to 
the  plate  bending  characteristics,  overall  equilibrium  is  closely 
sc  i .  n  t  a  i  ned  • 

The  composite  plate  element  can  be  used  in  the  existing 
-tutic  and  dynamic  analyses  contained  within  the  SAP  IV  program, 
i t  can  be  used  with  all  the  existing  elements  in  the  finite  ele¬ 
ment  library  as  long  as  the  model  effects  are  correct. 

The  element  was  not  being  developed  to  degenerate  to  a 
triangular  plate  element.  If  the  element  is  used  as  triangular, 
the  local  effect  is  "too  stiff".  If  the  fourth  of  the  quadri¬ 
lateral  nodes  is  placed  mid-plane  on  a  triangular  side,  a  better 
approximation  can  be  obtained. 

The  composite  clement  presently  does  not  contain  geometric 
tillering  elicits  ..uch  as  those  required  in  high  speed  centri- 
:  u,p.  machinery. 

further  Develoj .merit  s  : 

further  work  :  iiould  be  done  to  include  the  geometric  : 1 1  it  — 
inning  c  :  :  ic  i  cm.  ■  to  account  for  centripetal  acceleration 
..elects  of  high  spin  blade  systems.  This  geometric  matrix  of 
coefficients  would  allow  a  better  approximation  of  blade  bonding 
stresses  and  closer  representation  of  blade  natural  frequencies, 
at  high  spin. 

A  pro- processor  program  should  be  developed  to  handle  eomple 
material  ] ami nates  as  a  function  of  blade  position.  This  r  r<  - 
gram  would  compute  the  A,  B  and  D  matrices  and  thermal  vector: 
needed  in  the  SAP  program  and  produce  a  storage  file  for  stress 
recovery.  This  program  should  also  plot  all  information  i  or 
input  check  inf'. 


A  ..o-  should  Le  written  to  retrieve  stress  in- 

ion  divi  mater iai  information,  by  lamina,  to  be  used  with 
e: arm. tt ion  output  and  curvatures  at  the  mid-plane  surfaces  to 
.\-Juce  individual  lamina  stress  to  be  used  in  failure  criteria. 


i 
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appendices 


rEN.UX  A  -  liiPUT  TO  ELEMENT  TYPE  9  I:,  SA.-4A 

omplo  to  input  to  the  SAP  IV  prop, ram  ; «.  Ueucrio 
1.  T;.e  n.-w  element  TYPE  9  defined  in  f.A!  4 A  i 
e  .ini.i  mould  De  appended  to  reference  1.  The 
TYPE  9  i:-  similar  to  element  TYPE  6.  In  ;  act 
j  •  win.",  replace  TYPE;  f> .  The  format  o:  t.  e 
n  is  consistent  with  the  SAP  IV  input.  Varied, 
merer  must  be  right  justified  and  floating  no 
shoula  include  a  decimal. 


b'l 

J 


E  ELATE  INPUT  TO  SAP4A 


TYPE  9  Composite  Plate  Element  (QUADRILATERAL) 


Note  Columns  Variable 

A.  Control  Cara  (615,  110) 


5  NEAR ( 1 ) 

6-iC  NPAR ( 2 ) 

11-15  NPAR ( 3 ) 

(1)  16-20  NPAR ( 4 ) 


O  21-25  NPAR ( 5 ) 


2  6 -  3  C  NPAR ( 6 ) 

(3)  31-40  NPAR ( 7 ) 


Remark 


Number  9 

Number  of  Plate  Elements 
Number  of  different  materials 
Material  Type  Key 
=0  Composite  Material  Prop. 

Standard  Anisotropic  prop¬ 
erties  (same  as  TYPE  6) 

Number  of  Global  Material  Vectors 
If  zero  or  blank  then  global  X 
direction  is  assumed  to  be 
material  x  axis. 

Integration  Order  (default  set 
to  2 ) 

Rotation  Stiffness  Factor 
(integer  number) 


B.  Material  Property  Information 


Two  types  of  material  can  be 
input  to  element  type  9:  general 
composite  material  and  aniso¬ 
tropic  material. 

Composite  Material  Properties  ( NPAR ( 4 ) . EQ . 0 ) 

Five  cards  must  be  input  for 
every  different  material.  (NPAR(3)) 


Cara  1: 

(110, 

2  OX , 

4F10. 

0) 

j  -10 

NN 

Material  identification  number 

11-30 

Blank 

31-40 

DEN 

Mass  density 

41-50 

AT  ( 1 ) 

51-60 

AT  ( 2  ) 

Am 

thermal  vector 

component  s 

61-70 

AT  ( 3 ) 

—  1 

(equation  66) 

Card  2  : 

( 6F10 . 

0) 

1-10 

BT  ( 1 ) 

11-20 

BT  ( 2  ) 

21-30 

BT  (  3  ) 

-  i 

thermal  vector 

component : . 

31-40 

DT(  1) 

(equation  67) 

41-50 

DT(  2 ) 

I>T 

thermal  vector 

components 

51-60 

DT  (  3  ) 

—  i 

(equation  68) 

SR 


C.H-.1  3  : 


(6F10.0) 


:  - 1  o 

A  ( 1 , 1 ) 

11-2  0 

A(  1 , 2  ) 

A 

21-30 

A  ( 1 , 3  ) 

31-40 

A  ( 2 , 2  ) 

41-50 

A(  2  ,  3  ) 

51-60 

A(3  ,3) 

Card  4 : 

(6F10. 0) 

1-1C 

B  ( 1 , 1 ) 

11-20 

B(  1 , 2  ) 

B 

21-30 

B  ( 1 , 3  ) 

31-40 

B(2 ,2) 

41-50 

B(  2 , 3  ) 

51-60 

B  (  3 , 3  ) 

Card  5 : 

(6F10.0) 

1-10 

D  ( 1 , 1 ) 

11-20 

D( 1 , 2 ) 

D 

21-30 

D(1 , 3 ) 

31-40 

D(  2 , 2  ) 

41-50 

D(  2 , 3  ) 

51-60 

D(  3 , 3  ) 

Matrix  coefficients  (upper 
triangular)  (equation  B3) 


Matrix  coefficients  (upper 
triangular)  (equation  64) 


Matrix  coefficients  (upper 
triangular)  (equation  65) 


3.2  Anisotropic  Material  Properties  (NPAR( 4 ) . EQ . 1 ) 

(4)  Two  cards  must  be  input  for  every 

different  material  (NPAR(3)) 


Card  1 : 

(110,  2  OX , 

4F10.0) 

1-10 

NN 

Material  identification  number 

11-30 

Blank 

31-40 

DEN 

Mass  density 

41-50 

AX 

Thermal  expansion  coefficient  ax 

51-60 

AY 

Thermal  expansion  coefficient  ay 

61-70 

AXY 

Thermal  expansion  coefficient  axy 

Card  2  : 

(6F10.0) 

1-10 

CXX 

Elasticity  element  Cxx 

11-20 

CXY 

Elasticity  element  Cxy 

21-30 

CXS 

Elasticity  element  Cxs 

31-40 

CYY 

Elasticity  element  Cyy 

41-50 

CYS 

Elasticity  element  Cys 

51-60 

GXY 

Elasticity  element  Gxy 

59 


Material 

Vectors  (25, 

5X ,  3F10.0) 

NPAR(S) 

material  ve 

ctors  must  be 

input  (except  if 

zero  or  blank) 

1-5 

NV 

Material  vector  identification 

number 

6-1C 

Blank 

11-20 

DX 

X  direction 

cosine 

21-30 

DY 

Y  direction 

cosine 

31-40 

DZ 

Z  direction 

cosine 

Element  Load  Multipliers  (5  cards) 
Card  1:  (4F10.0) 


1-10 

PA 

Distributed  lateral  load  multi- 
plier  for  load  case  A 

11-20 

P5 

Distributed  lateral  load  multi¬ 
plier  for  load  case  B 

21-30 

PC 

Distributed  lateral  load  multi¬ 
plier  for  load  case  C 

3 1  —  “f  C 

i'D 

Distributed  lateral  load  multi¬ 
plier  for  load  case  D 

Card  2  : 

(4F10.0) 

1-10 

TA 

Temperature  multiplier  for  load 
case  A 

11-20 

TB 

Temperature  multiplier  for  load 

case  B 

21-30 

TC 

Temperature  multiplier  for  load 
case  C 

31-40 

TD 

Temperature  multiplier  for  load 
case  D 

Card  3 : 

(4F10.0) 

1-10 

XA 

X-direction  acceleration  for 
load  case  A 

11-20 

XB 

X-direction  acceleration  for 
load  case  B 

21-30 

XC 

X-direction  acceleration  for 
load  case  C 

31-40 

XD 

X-direction  acceleration  for 
load  case  D 

60 


Card  4 : 

( * i  FI  0.0) 

1-10 

V  A 

Y-di rect ion 
load  case  A 

aeee  1  .-rat  ior. 

!V,;- 

1 

t'O 

o 

Yn 

Y-n i rection 

accelerat ion 

r  c>:  ■ 

iOtid  p 

21-30 

YC 

Y-di rect i on 

acce  1  or.,  t  ior. 

!  1 1 1 ' 

load  e.ir.i-  C 

31-10 

V  • 

Y-d  i  rect  1  or; 
load  case  D 

acce  r e rat  ion 

:  cr 

Card  S  : 

(4F1G. 0) 

1-10 

3  A 

Z-dxrection 
load  case  A 

acce lerat ion 

for 

11-2  0 

7,3 

Z-direction 
load  case  E 

acce 1  oration 

fc.r 

tv' 

l—1 

1 

C*J 

o 

ZC 

Z-direction 
load  case  C 

acceleration 

:  o  r 

31-4  0 

ZD 

Z-direction 

acceleration 

f  or 

load  case  D 


s .  1 

lament  Cards  ( 

515,  12, 

13,  12,  13,  15,  4F10.0) 

One  card 

for  each 

h  i  'AR  (  2  )  element. 

1-5 

NN 

Element  number 

(5) 

6-  i  0 

I 

Node  I 

11-  lb 

J 

•Node  1 1 

L  '1 

i.  O  —  z. 

K 

Node  K 

2  i  -  ?  S 

i. 

Node  I, 

.  ;  ) 

2  G  -  2  7 

NO 

ho.  or  Gauss  integration  ; o i n : . 

(7) 

2  8  -  3  3 

IV 

Material  vector  identification 
number 

(d) 

3  I  -  3  2 

iRhdSt: 

Previous  Element  re-uae  coco 

-0  new  element 

-1  use  previous  element 

(9) 

3  3- 3b 

T,M 

Material  identification  number 

(13) 

36-40 

I NCL 

Element  generation  parameter 

(11  ) 

41-50 

TH 

Element  thickness 

(12) 

51-60 

PR 

Element  lateral  pressure 

(13) 

61-70 

TO 

Mean  temperature  variation  from 
the  reference  level  in  uncle  - 
! ormed  pos i t i on . 

71-80 

TO 

Moan  temperature  gradient  acres 

t  ho  she  11  till  ekness . 

h< .  to;. : 

(1)  clement  TYI'i,  4  allow:;  two  *  i  i  t  t  eren  I  i  o.-m  input.  ’i'he 

:  irst  form  is  i  or  laminate  matrices  whi  »<•  t  ho  „eeond  is  1  <>t 
anisotropic  matrices .  The  Later  form  is  ilcnric.il  to  ele¬ 
ment  TYPE  6  input. 


hi 


(2)  A  nuiU-Ti.il  plobal  axis  must  be  defined  relative  to  the 
material  properties  formed  in  the  A,  B  and  D  matrices. 

(3)  Rotational  Stillness  Factor  is  set  by  multiplying,  NI'AR(V) 
times  l.F-8.  default  for’  NFAK(7)  is  100. 

(4)  Material  input  in  this  section  is  identical  to  that  of 
element  TYPL  6. 

(3)  The  T,J,K  and  I,  indices  del  ine  tlie  element  connectivity 
and  also  the  element  normal.  The  element  "z"  coordinate 
is  !  or'med  bv  the  ripht  hand  rule  as  T  peer,  to  d  p.oes  to  K, 
etc.  The  element  local  axis  is  determined  by  the  protec¬ 
tion  of  the  plobal  material  axis  onto  the  element.  Once 
the  x-axls  is  determined,  the  local  y  is  formed  from  z 
and  x.  All  stress  output  is  in  this  re  I o pence.  If  node  L 
espial  K  or  i!  zero  or  left  blank,  the  propram  will  assume 
that  the  element  is  trianpulur.  The  resultinp  local  stiff¬ 
ness  is  then  "too  stiff". 

( >  )  The  number  ot  Gauss  intep.ration  points  can  vary  as  2  or  3. 

It  a  value  is  set  above  or  below  these  numbers,  the  prop ram 
will  reset  it  to  ?.  The  del  au  1 1  value?  is-  set  ft  PAR  ( 6 ) . 

(7)  The  plobal  material  vector-  must  be  preater  than  or  equal 

to  1  and  less,  than  or  equal  to  N  PAR  ( 5 ) .  If  N1’AR(5)  is  blank 
•  >r  zero,  then  NI’AR(  5 )  is  set  e<)ua  l  to  1  and  the  plobal  vec¬ 
tor  is.  al  i  piled  alonp.  the  plobal  X  axis..  be  fault  value  is 
set  to  1 . 

(3)  if  an  element  has.  the  same  planar  size,  same  oriental  ic-n  in 
space  and  the  same  element  load  i  np  parameters  as  the  pre¬ 
vious.  f-hwnl  ,  then  sett  inn  1  KFUGh  equal  to  1  will  use  the 
same  plobal  element  stiffness  and  load  vectors  for  assembly, 
default  if’-  set  to  0. 

(3)  The  material  Id  number  must  be  between  1  and  NFAR(3).  de¬ 
fault  is  set  to  1. 

(Ill)  element  Generation  barometer:  Moment  cards  must  be  in 

. •  1  emerit  number  sequence.  ll  clement  cards-  are  omitted,  the 
propr-am  wi  I  1  pi-nf-r.i  te  the  mis.sinp  ear-ds.  ,is  follows: 

Tli*  •  increment  I  or  the-  element  number  i  one. 

i  .  e .  ,  UM .  ,  -  NfJ  •  +  1 

Tli.  •  r  ■( , t r*  ■ 1  |  ■<  'll- 1  i  u i-  i  nc  r<  -m.  ■  n  t  I  >  > r  nodal  <  i  .n  nee  I  i  v  i  I  v  i  s  1  N' b  I . . 


•  > 


M  +  l 


Ii  +  INCL 


J.+i  = 

J.  + 

l 

INCL 

Ki+i  = 

K.  + 

l 

INCL 

Li+1  = 

L.  + 

l 

INCL 

If  INCL  is  left  blank 

then 

INCL  is 

set  to  1.  Material 

identification,  elemen 

t  thi 

ckness , 

distributed  lateral 

load,  temperature  and  temperature  gradient  for  the  element 
are  then  the  same  as  for  the  first  element  in  the  generated 
group.  The  last  element  card  must  be  input  to  exit  element 
group  properly. 

(11)  The  plate  thickness  is  used  mainly  to  compute  the  mass  of 
the  element.  Default  is  set  to  1.0. 

(12)  The  pressure  is  normal  to  the  surface  of  the  element.  The 
positive  direction  for  the  pressure  loading  vector  is  in 
the  positive  direction  of  the  local  z  coordinate. 

(13)  The  temperature  required  is  the  mean  temperature  difference 
(TO)  from  the  reference  temperature  of  the  element  in  a 
undeformed  state.  The  TG  is  the  mean  thermal  gradient 
through  the  element  thickness. 
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APPENDIX  B 


INPUT  TO  PRE-PROCESSOR  PROGRAM  LAYUP' 


A  pre-proct  ssor  program  wa:.  developed  to  calculate  the 
and  D  natrice.i  ts  described  in  equations  63  through  63  and 
t  e  thermal  load  vec  i  ors  describee  in  equation:-  66  throupn  68 . 
'.tie  material  Lnt  ornat ion  by  layer  in  input  no  lamina  position 
and  : iber  orientation.  Tne  pre-processor  then  computes  the  0 
it:  it  rices  invoivin,  the  tensor  transformation  and  constructs  the 
element  material  matrix  E  relative  to  the  mid-plane  of  the 
plate.  Similar  thermal  vectors  are  also  computed  and  output. 
This  information  is  then  used  directly  for  SAP4A  -  TYPE  9. 
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INPUT  TO  PROGRAM  LAYUP 


A. 

Number  of  Cases 

Card  (12) 

j  Columns 

Variable 

Remark 

CD 

1-2 

N CASES 

Enter  the  total  number  of 
laminate  configurations  to  be 
considered . 

5 . 

Heading  or  Titl 

e  Card  (8A10) 

(2) 

1-60 

ITITLE 

Enter  the  title  information  to 
be  printed  with  the  output. 

C. 

Laminate  Data  Card  (415) 

(3) 

1-0 

NLAM 

Enter  the  number  of  laminae  in 
the  laminate. 

(‘  t) 

6-IC 

NMAT 

Enter  the  number  of  different 
materials  in  the  laminate. 

11-15 

ITRAN 

If  ITRAN  is  zero  or  blank  the 
transformed  thermal  properties 
of  each  lamina  will  not  be  part 
of  the  output. 

16-20 

I0RD 

If  IORD  is  left  blank  the  lami¬ 
nate  ordinates  will  not  be  output 

0 . 

Material  Proper 

ty  Card  (6F10. 

3) 

(b) 

1-10 

Ell(K) 

Enter  lamina  modulus  in  fiber 
direction 

"  1  —  >  ,'j 

L2  2 ( K ) 

Enter  lamina  modulus,  in  direction 
transverse  to  fibers. 

: :  -  j  o 

XN1(K) 

Enter  maior  Poisson's  ratio. 

3  - -40 

ni2(K) 

Enter  lamina  shear  modulus . 

l>  1  -50 

THERM 1 ( K) 

Enter  C.T.l.  in  : iber  direction. 

51-60 

THERM 2 ( K) 

Enter  C.T.E.  in  transverse 
direction . 

L . 

earn in a  Data  Car 

d  (no. 3,  no. 

F10. 3) 

i- 1  e 

T  ( J ) 

Enter  thickness  of  the  lamina 

11-20 

MATL(J) 

Enter  the  number  that t^dent if ies 
the  material  of  the  E  ‘  lamina^w 
Enter  the  orientation  of  the  J 
lamina  with  respect  to  the 
laminate  axes. 

2  1-30 

PJI(J) 

Notes : 

(1)  There  are  no  program  restrictions  on  the  number  of  cases 
that  may  be  analyzed  in  a  single  run. 

(2)  Begin  each  new  data  case  with  a  heading  card. 

(3)  The  program  is  currently  capable  of  handling  up  to  48  lami¬ 
nae  per  layup.  This  can  be  increased  by  changing  the 
appropriate  dimension  statements  as  shown  in  the  program 
listing. 

(4)  Although  no  limitation  on  the  number  of  different  materials 
need  be  imposed,  the  program  dimension  statements  currently 
allow  for  a  maximum  of  NMAT  =  4.  However,  this  can  also  be 
increased  if  necessary. 

(5)  This  material  data  is  vendor  information.  One  card  is 
required  for  each  material  (i.e.,  K  =  1,  NMAT). 

(6)  One  card  is  required,  per  layer,  in  the  laminate.  (i.e., 

J  =  1,  NLAM) . 


APPENDIX  C 


COMPUTER  RUN 


INPUT  AND  OUTPUT 


This  appendix  contains  the  complete  inpul  and  resulting 
output  of  the  SAP4A  program  using  the  new  composite  plate  ele¬ 
ment  (TYPE  9).  The  example  is  the  model  contained  in  section 
4.7  (Curved  Plate  under  Uniform  Pressure).  The  first  section 
contains  a  listing  of  the  card  images  used  to  execute  the  pro¬ 
gram.  The  second  section  is  the  complete  output  resulting 
from  this  input. 
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